Interfacing with Lattice Problems, Integrators, and Solutions

We have previously described how to retrieve species and parameter values stored in non-spatial problems, integrators, and solutions. This section describes similar workflows for simulations based on LatticeReactionSystems.

Generally, while for non-spatial systems these operations can typically be done by indexing a structure directly, e.g. through

sol[:X]

there are no equally straightforward interfaces for spatial simulations. Typically, helper functions have to be used, e.g

lat_getu(sol, :X, lrs)

Furthermore, there are some cases of interfacing which are currently not supported (e.g. updating parameter values in JumpProblems). It is likely that these interfaces will be improved in the future (i.e. by introducing a similar syntax to the current non-spatial one). Finally, we note that many of the functions presented below have not been as extensively optimised for performance as other parts of the Catalyst package. Hence, you should take care when designing workflows which requires using them a large number of times.

Retrieving values from lattice simulations

Let us consider a simulation of a LatticeReactionSystem:

using Catalyst, OrdinaryDiffEqDefault
two_state_model = @reaction_network begin
    (k1,k2), X1 <--> X2
end
diffusion_rx = @transport_reaction D X1
lattice = CartesianGrid((2,3))
lrs = LatticeReactionSystem(two_state_model, [diffusion_rx], lattice)

u0 = [:X1 => [0.0 0.0 0.0; 2.0 2.0 2.0], :X2 => 0.0]
tspan = (0.0, 1.0)
ps = [:k1 => 2.0, :k2 => 1.0, :D => 0.1]
oprob = ODEProblem(lrs, u0, tspan, ps)
sol = solve(oprob)

To retrieve the values of $X1$ across the simulation we use the lat_getu function. It takes three arguments:

• The solution objects from which we wish to retrieve values.
• The species which values we wish to retrieve.
• The LatticeReactionSystem which was simulated.

lat_getu(sol, :X1, lrs)

13-element Vector{Matrix{Float64}}:
 [0.0 0.0 0.0; 2.0 2.0 2.0]
 [4.988650432250024e-5 4.988650432250024e-5 4.988650432250024e-5; 1.998952234018206 1.998952234018206 1.998952234018206]
 [0.0005458858864143923 0.0005458858864143923 0.0005458858864143925; 1.9885184230183552 1.9885184230183552 1.9885184230183552]
 [0.004183579842347592 0.004183579842347592 0.004183579842347594; 1.9110504558577903 1.9110504558577903 1.9110504558577903]
 [0.01055512495082666 0.01055512495082666 0.010555124950826672; 1.7708833918940514 1.7708833918940514 1.7708833918940514]
 [0.0175113054328957 0.017511305432895696 0.017511305432895714; 1.6099324445575665 1.6099324445575665 1.6099324445575667]
 [0.02475470074637461 0.02475470074637461 0.02475470074637463; 1.430770723996124 1.4307707239961238 1.430770723996124]
 [0.031518542470780904 0.03151854247078091 0.031518542470780925; 1.2486002352246621 1.248600235224662 1.2486002352246621]
 [0.03729771788004629 0.03729771788004631 0.0372977178800463; 1.0766135558411272 1.0766135558411267 1.0766135558411272]
 [0.04182477924153854 0.04182477924153853 0.04182477924153854; 0.9284335024464289 0.9284335024464287 0.9284335024464289]
 [0.045285465842461646 0.045285465842461625 0.045285465842461646; 0.8115561309730012 0.8115561309730013 0.8115561309730012]
 [0.04818220881966842 0.048182208819668414 0.04818220881966845; 0.7274838954502673 0.7274838954502675 0.7274838954502671]
 [0.05042586336228593 0.05042586336228592 0.05042586336228594; 0.6826244818846041 0.6826244818846042 0.6826244818846041]

Here, the output is a vector with $X1$'s value at each simulation time step. How the species's value is represented at each time step depends on the lattice which was originally used to create the LatticeReactionSystem:

• For Cartesian lattices, an array of the same size as the Cartesian lattice is used. Each array element corresponds to the species's value in the corresponding compartment.
• For masked lattices, a sparse array of the same size as the masked lattice is used. Each filled array element corresponds to the species's value in the corresponding compartment. Unfilled array elements correspond to positions without compartments.
• For unstructured (graph) lattices, vectors are used. The i'th element in the vectors corresponds to the species's value in the i'th compartment.

Unlike for non-spatial simulations, lat_getu does not take vector (e.g. lat_getu(sol, [:X1, :X2], lrs)) or symbolic expression (e.g. lat_getu(sol, [X1 + X2], lrs)) inputs. However, it is possible to use symbolic variables as input (e.g. lat_getu(sol, two_state_model.X1, lrs)).

Retrieving lattice simulations values at specific time points

Just like for non-spatial solutions, it is possible to access the simulation's values at designated time points. This is possible even if the simulation did not stop at those specific time points (in which case an interpolated value is returned). To do this, the desired time points to sample are provided as a vector to lat_getu using the optional argument t. E.g. here we retrieve the simulation's (interpolated) values at time points 0.5 and 0.75:

lat_getu(sol, :X1, lrs; t = [0.5, 0.75])

2-element Vector{Matrix{Float64}}:
 [0.04200952741591302 0.04200952741591301 0.04200952741591302; 0.9221641543216159 0.9221641543216158 0.9221641543216159]
 [0.04694938928015799 0.046949389280157974 0.04694938928015798; 0.760247120727087 0.760247120727087 0.760247120727087]

Retrieving and updating species values in problems and integrators

Let us consider a spatial ODEProblem

using Catalyst, OrdinaryDiffEqDefault
two_state_model = @reaction_network begin
    (k1,k2), X1 <--> X2
end
diffusion_rx = @transport_reaction D X1
lattice = CartesianGrid((2,3))
lrs = LatticeReactionSystem(two_state_model, [diffusion_rx], lattice)

u0 = [:X1 => [0.0 0.0 0.0; 2.0 2.0 2.0], :X2 => 0.0]
tspan = (0.0, 1.0)
ps = [:k1 => 2.0, :k2 => 1.0, :D => 0.1]
oprob = ODEProblem(lrs, u0, tspan, ps)

We can retrieve the species values stored in oprob using the lat_getu function. It uses identical syntax as for simulations (except that you cannot specify a time point). However, it returns a single set of species values (while for simulations it returns a vector across different time steps):

lat_getu(oprob, :X1, lrs)

2×3 Matrix{Float64}:
 0.0  0.0  0.0
 2.0  2.0  2.0

Again, the format used corresponds to the lattice used to create the original LatticeReactionSystem. Here, even if a species has homogeneous values, the full format is used.

lat_getu(oprob, :X2, lrs)

2×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0

For both problems and integrators, species values can be updated using the lat_setu! function. It uses a similar syntax as lat_getu, but takes a fourth argument which is the new values to use for the designated species:

lat_setu!(oprob, :X1, lrs, [1.0 2.0 3.0; 4.0 5.0 6.0])

Here, the same format (which depends on the used lattice) is used for the species's new values, as which is used when initially designating their initial conditions. I.e. to make $X1$'s initial condition values uniform we can call

lat_setu!(oprob, :X1, lrs, 1.0)

Species values in integrators can be interfaced with using identical syntax as for problems.

Retrieving and updating parameter values in problems and integrators

Retrieval and updating of parameter values for problems and integrators works similarly as for species, but with the following differences:

• The lat_getp and lat_setp! functions are used.
• It is currently not possible to interface with parameter values of JumpProblems and their integrators.
• After parameter values are modified, the rebuild_lat_internals! function must be applied before the problem/integrator can be used for further analysis.
• Updating of edge parameters is limited and uses a different interface.

Let us consider the spatial ODEProblem we previously declared. We can check the value of $k1$ by using lat_getp

lat_getp(oprob, :k1, lrs)

2×3 Matrix{Float64}:
 2.0  2.0  2.0
 2.0  2.0  2.0

Next, we can update it using lat_setp! (here we also confirm that it now has the updated values):

lat_setp!(oprob, :k1, lrs, [1.0 2.0 3.0; 4.0 5.0 6.0])
lat_getp(oprob, :k1, lrs)

2×3 Matrix{Float64}:
 1.0  2.0  3.0
 4.0  5.0  6.0

If we now were to simulate oprob, the simulation would not take the updated value of $k1$ into account. For our changes to take effect we might first need to call rebuild_lat_internals! with oprob as an input

rebuild_lat_internals!(oprob)

There are two different circumstances when rebuild_lat_internals! must be called:

• When modifying the value of an edge parameter.
• When changing a parameter from having spatially uniform values to spatially non-uniform values, or the other way around.

Parameter values of integrators can be interfaced with just like for problems (this is primarily relevant when using callbacks). Again, after doing so, the rebuild_lat_internals! function might need to be applied to the integrator.

Retrieving and updatingedge parameter values in problems and integrators

The lat_getp and lat_setp! functions cannot currently be applied to edge parameters. Instead, to access the value of an edge parameter, use

oprob.ps[:D]

1×1 SparseArrays.SparseMatrixCSC{Float64, Int64} with 1 stored entry:
 0.1

To update an edge parameter's value, use

oprob.ps[:D] = [0.2]

This interface is somewhat limited, and the following aspects should be noted:

• Edge parameter values can only be interfaced with if the edge parameter's value is spatially uniform.
• When accessing an (spatially uniform) edge parameter's value, its single value will be encapsulated in a vector.
• When setting an (spatially uniform) edge parameter's value, you must encapsulate the new value in a vector.